GLM Example 1

Author

Murray Logan

Published

30/07/2023

1 Preparations

Load the necessary libraries

library(glmmTMB)       #for model fitting
library(car)           #for regression diagnostics
library(broom)         #for tidy output
library(ggfortify)     #for model diagnostics
library(DHARMa)        #for residual diagnostics
library(see)           #for plotting residuals
library(knitr)         #for kable
library(effects)       #for partial effects plots
library(ggeffects)     #for partial effects plots
library(emmeans)       #for estimating marginal means
library(modelr)        #for auxillary modelling functions
library(tidyverse)     #for data wrangling
library(lindia)        #for diagnostics of lm and glm
library(performance)   #for residuals diagnostics
library(sjPlot)        #for outputs
library(report)        #for reporting methods/results
library(easystats)     #framework for stats, modelling and visualisation

2 Scenario

Here is an example from Fowler, Cohen, and Jarvis (1998). An agriculturalist was interested in the effects of fertiliser load on the yield of grass. Grass seed was sown uniformly over an area and different quantities of commercial fertiliser were applied to each of ten 1 m² randomly located plots. Two months later the grass from each plot was harvested, dried and weighed. The data are in the file fertilizer.csv in the data folder.

Figure 1: Field of grass
Table 1: Format of the fertilizer.csv data file
FERTILIZER YIELD
25 84
50 80
75 90
100 154
125 148
... ...
Table 2: Description of the variables in the fertilizer data file
FERTILIZER: Mass of fertiliser (g/m²) - Predictor variable
YIELD: Yield of grass (g/m²) - Response variable

The aim of the analysis is to investigate the relationship between fertiliser concentration and grass yield.

3 Read in the data

We will start off by reading in the Fertiliser data. There are many functions in R that can read in a CSV file. We will use a the read_csv() function as it is part of the tidyverse ecosystem.

fert <- read_csv("../public/data/fertilizer.csv", trim_ws = TRUE)
Rows: 10 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (2): FERTILIZER, YIELD

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

After reading in a dataset, it is always a good idea to quickly explore a few summaries in order to ascertain whether the imported data are correctly transcribed. In particular, we should pay attention to whether there are any unexpected missing values and ensure that each variable (column) has the expected class (e.g. that variables we expected to be considered numbers are indeed listed as either <dbl> or <int> and not <char>).

glimpse(fert)
Rows: 10
Columns: 2
$ FERTILIZER <dbl> 25, 50, 75, 100, 125, 150, 175, 200, 225, 250
$ YIELD      <dbl> 84, 80, 90, 154, 148, 169, 206, 244, 212, 248
## Explore the first 6 rows of the data
head(fert)
# A tibble: 6 × 2
  FERTILIZER YIELD
       <dbl> <dbl>
1         25    84
2         50    80
3         75    90
4        100   154
5        125   148
6        150   169
str(fert)
spc_tbl_ [10 × 2] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
 $ FERTILIZER: num [1:10] 25 50 75 100 125 150 175 200 225 250
 $ YIELD     : num [1:10] 84 80 90 154 148 169 206 244 212 248
 - attr(*, "spec")=
  .. cols(
  ..   FERTILIZER = col_double(),
  ..   YIELD = col_double()
  .. )
 - attr(*, "problems")=<externalptr> 
fert |> datawizard::data_codebook()
fert (10 rows and 2 variables, 2 shown)

ID | Name       | Type    | Missings |    Values |  N
---+------------+---------+----------+-----------+---
1  | FERTILIZER | numeric | 0 (0.0%) | [25, 250] | 10
---+------------+---------+----------+-----------+---
2  | YIELD      | numeric | 0 (0.0%) | [80, 248] | 10
-----------------------------------------------------

4 Exploratory data analysis

The individual responses (\(y_i\), observed yields) are each expected to have been independently drawn from normal (Gaussian) distributions (\(\mathcal{N}\)). These distributions represent all the possible yields we could have obtained at the specific (\(i^{th}\)) level of Fertiliser. Hence the \(i^{th}\) yield observation is expected to have been drawn from a normal distribution with a mean of \(\mu_i\).

Although each distribution is expected to come from populations that differ in their means, we assume that all of these distributions have the same variance (\(\sigma^2\)).

Model formula: \[ \begin{align} y_i &\sim{} \mathcal{N}(\mu_i, \sigma^2)\\ \mu_i &= \beta_0 + \beta_1 x_i \end{align} \]

where \(y_i\) represents the \(i\) observed values, \(\beta_0\) and \(\beta_1\) represent the intercept and slope respectively, and \(\sigma^2\) represents the estimated variance.

We intend to explore the relationship between Yield and Fertiliser using simple linear regression. Such an analysis assumes:

  • Normality: each population is assumed to be normally distributed. This can be initially assessed by exploring the distribution of the response via either boxplots, violin plots or histograms (depending on sample size).
  • Homogeneity of variance: each population is assumed to be equally varied. This can be initially assessed by exploring the distribution of observed values around an imaginary line of best fit through the cloud of data formed by a scatterplot of Yield against Fertiliser. If the spread of points increases (or decreases) along the trend line, it is likely that the populations are not equally varied.
  • Linearity: in fitting a straight line through the cloud of data, we are assuming that a linear trend is appropriate. If the trend is not linear, then the inferred relationship may be miss-representative of the true relationship. In addition to an imaginary line, we could fit either a loess or linear smoother through the cloud to help us assess the likelihood of non-linearity.
  • Independence: to help ensure that the estimates are unbiased, we must assume that all observations are independent. In this case, we assume that the observations were collected from a randomised design in which the level of Fertiliser was applied randomly to the different patches of the field. If however, the researchers had simply moved along the field applying progressively higher Fertiliser concentrations, then there is a chance that they have introduced additional biases or confounding factors. Similarly, if the Fertiliser levels were applied in increasing doses over time, there might also be biases. In the absence of any spatial or temporal information associated with the collection of data, we cannot directly assess this assumption.
ggplot(fert, aes(y = YIELD, x = FERTILIZER)) +
  geom_point() +
  geom_smooth()

ggplot(fert, aes(y = YIELD, x = FERTILIZER)) +
  geom_point() +
  geom_smooth(method = "lm")

ggplot(fert, aes(y = YIELD)) +
    geom_boxplot(aes(x = 1))

ggplot(fert, aes(y = YIELD)) +
  geom_violin(aes(x = 1))

ggplot(fert, aes(x = YIELD)) +
  geom_histogram()

ggplot(fert, aes(x = YIELD)) +
  geom_density()

Conclusions:

  • there is no evidence of non-normality, non-homogeneity of variance or non-linearity
  • we have no initial reason to suspect that the assumptions will not be satisfied and thus it is reasonable to assume that the results will be reliable.

5 Fit the model

fert_lm <- lm(YIELD ~ 1 + FERTILIZER, data = fert)
fert_lm <- lm(YIELD ~ FERTILIZER, data = fert)

The lm() function returns a list comprising the estimated coefficients, the residuals, the data used to fit the model and various other attributes associated with the model fitting procedure.

## Get the name of the attributes within the lm model
attributes(fert_lm)
$names
 [1] "coefficients"  "residuals"     "effects"       "rank"         
 [5] "fitted.values" "assign"        "qr"            "df.residual"  
 [9] "xlevels"       "call"          "terms"         "model"        

$class
[1] "lm"
## Explore all the conttents of the fitted model
str(fert_lm)
List of 12
 $ coefficients : Named num [1:2] 51.933 0.811
  ..- attr(*, "names")= chr [1:2] "(Intercept)" "FERTILIZER"
 $ residuals    : Named num [1:10] 11.78 -12.5 -22.79 20.93 -5.36 ...
  ..- attr(*, "names")= chr [1:10] "1" "2" "3" "4" ...
 $ effects      : Named num [1:10] -517.03 184.25 -23.73 18.66 -8.96 ...
  ..- attr(*, "names")= chr [1:10] "(Intercept)" "FERTILIZER" "" "" ...
 $ rank         : int 2
 $ fitted.values: Named num [1:10] 72.2 92.5 112.8 133.1 153.4 ...
  ..- attr(*, "names")= chr [1:10] "1" "2" "3" "4" ...
 $ assign       : int [1:2] 0 1
 $ qr           :List of 5
  ..$ qr   : num [1:10, 1:2] -3.162 0.316 0.316 0.316 0.316 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:10] "1" "2" "3" "4" ...
  .. .. ..$ : chr [1:2] "(Intercept)" "FERTILIZER"
  .. ..- attr(*, "assign")= int [1:2] 0 1
  ..$ qraux: num [1:2] 1.32 1.27
  ..$ pivot: int [1:2] 1 2
  ..$ tol  : num 1e-07
  ..$ rank : int 2
  ..- attr(*, "class")= chr "qr"
 $ df.residual  : int 8
 $ xlevels      : Named list()
 $ call         : language lm(formula = YIELD ~ FERTILIZER, data = fert)
 $ terms        :Classes 'terms', 'formula'  language YIELD ~ FERTILIZER
  .. ..- attr(*, "variables")= language list(YIELD, FERTILIZER)
  .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:2] "YIELD" "FERTILIZER"
  .. .. .. ..$ : chr "FERTILIZER"
  .. ..- attr(*, "term.labels")= chr "FERTILIZER"
  .. ..- attr(*, "order")= int 1
  .. ..- attr(*, "intercept")= int 1
  .. ..- attr(*, "response")= int 1
  .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. ..- attr(*, "predvars")= language list(YIELD, FERTILIZER)
  .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
  .. .. ..- attr(*, "names")= chr [1:2] "YIELD" "FERTILIZER"
 $ model        :'data.frame':  10 obs. of  2 variables:
  ..$ YIELD     : num [1:10] 84 80 90 154 148 169 206 244 212 248
  ..$ FERTILIZER: num [1:10] 25 50 75 100 125 150 175 200 225 250
  ..- attr(*, "terms")=Classes 'terms', 'formula'  language YIELD ~ FERTILIZER
  .. .. ..- attr(*, "variables")= language list(YIELD, FERTILIZER)
  .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:2] "YIELD" "FERTILIZER"
  .. .. .. .. ..$ : chr "FERTILIZER"
  .. .. ..- attr(*, "term.labels")= chr "FERTILIZER"
  .. .. ..- attr(*, "order")= int 1
  .. .. ..- attr(*, "intercept")= int 1
  .. .. ..- attr(*, "response")= int 1
  .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. .. ..- attr(*, "predvars")= language list(YIELD, FERTILIZER)
  .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
  .. .. .. ..- attr(*, "names")= chr [1:2] "YIELD" "FERTILIZER"
 - attr(*, "class")= chr "lm"
## Return the data used to fit the model 
fert_lm$model
   YIELD FERTILIZER
1     84         25
2     80         50
3     90         75
4    154        100
5    148        125
6    169        150
7    206        175
8    244        200
9    212        225
10   248        250
## Return the estimated coefficients
fert_lm$coefficients
(Intercept)  FERTILIZER 
 51.9333333   0.8113939 

Rather than directly access these attributes (which may have different names in different modelling functions), there are extraction functions that reach into the model output and extract the information in a convenient way.

## Extract the estimated coefficients
coef(fert_lm)
(Intercept)  FERTILIZER 
 51.9333333   0.8113939 
## Extract the fitted (predicted) values
fitted(fert_lm)
        1         2         3         4         5         6         7         8 
 72.21818  92.50303 112.78788 133.07273 153.35758 173.64242 193.92727 214.21212 
        9        10 
234.49697 254.78182 
## Extract the residules
resid(fert_lm)
         1          2          3          4          5          6          7 
 11.781818 -12.503030 -22.787879  20.927273  -5.357576  -4.642424  12.072727 
         8          9         10 
 29.787879 -22.496970  -6.781818 
args(residuals.lm)
function (object, type = c("working", "response", "deviance", 
    "pearson", "partial"), ...) 
NULL
fert_lm1 <- lm(YIELD ~ center(FERTILIZER), data = fert)

The lm() function returns a list comprising the estimated coefficients, the residuals, the data used to fit the model and various other attributes associated with the model fitting procedure.

## Get the name of the attributes within the lm model
attributes(fert_lm1)
$names
 [1] "coefficients"  "residuals"     "effects"       "rank"         
 [5] "fitted.values" "assign"        "qr"            "df.residual"  
 [9] "xlevels"       "call"          "terms"         "model"        

$class
[1] "lm"
## Explore all the conttents of the fitted model
str(fert_lm1)
List of 12
 $ coefficients : Named num [1:2] 163.5 0.811
  ..- attr(*, "names")= chr [1:2] "(Intercept)" "center(FERTILIZER)"
 $ residuals    : Named num [1:10] 11.78 -12.5 -22.79 20.93 -5.36 ...
  ..- attr(*, "names")= chr [1:10] "1" "2" "3" "4" ...
 $ effects      : Named num [1:10] -517.03 184.25 -23.73 18.66 -8.96 ...
  ..- attr(*, "names")= chr [1:10] "(Intercept)" "center(FERTILIZER)" "" "" ...
 $ rank         : int 2
 $ fitted.values: Named num [1:10] 72.2 92.5 112.8 133.1 153.4 ...
  ..- attr(*, "names")= chr [1:10] "1" "2" "3" "4" ...
 $ assign       : int [1:2] 0 1
 $ qr           :List of 5
  ..$ qr   : num [1:10, 1:2] -3.162 0.316 0.316 0.316 0.316 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:10] "1" "2" "3" "4" ...
  .. .. ..$ : chr [1:2] "(Intercept)" "center(FERTILIZER)"
  .. ..- attr(*, "assign")= int [1:2] 0 1
  ..$ qraux: num [1:2] 1.32 1.27
  ..$ pivot: int [1:2] 1 2
  ..$ tol  : num 1e-07
  ..$ rank : int 2
  ..- attr(*, "class")= chr "qr"
 $ df.residual  : int 8
 $ xlevels      : Named list()
 $ call         : language lm(formula = YIELD ~ center(FERTILIZER), data = fert)
 $ terms        :Classes 'terms', 'formula'  language YIELD ~ center(FERTILIZER)
  .. ..- attr(*, "variables")= language list(YIELD, center(FERTILIZER))
  .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. ..$ : chr [1:2] "YIELD" "center(FERTILIZER)"
  .. .. .. ..$ : chr "center(FERTILIZER)"
  .. ..- attr(*, "term.labels")= chr "center(FERTILIZER)"
  .. ..- attr(*, "order")= int 1
  .. ..- attr(*, "intercept")= int 1
  .. ..- attr(*, "response")= int 1
  .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. ..- attr(*, "predvars")= language list(YIELD, center(FERTILIZER, center = 137.5, verbose = FALSE))
  .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
  .. .. ..- attr(*, "names")= chr [1:2] "YIELD" "center(FERTILIZER)"
 $ model        :'data.frame':  10 obs. of  2 variables:
  ..$ YIELD             : num [1:10] 84 80 90 154 148 169 206 244 212 248
  ..$ center(FERTILIZER): 'dw_transformer' num [1:10] -112.5 -87.5 -62.5 -37.5 -12.5 ...
  .. ..- attr(*, "center")= num 138
  .. ..- attr(*, "scale")= num 1
  .. ..- attr(*, "robust")= logi FALSE
  ..- attr(*, "terms")=Classes 'terms', 'formula'  language YIELD ~ center(FERTILIZER)
  .. .. ..- attr(*, "variables")= language list(YIELD, center(FERTILIZER))
  .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:2] "YIELD" "center(FERTILIZER)"
  .. .. .. .. ..$ : chr "center(FERTILIZER)"
  .. .. ..- attr(*, "term.labels")= chr "center(FERTILIZER)"
  .. .. ..- attr(*, "order")= int 1
  .. .. ..- attr(*, "intercept")= int 1
  .. .. ..- attr(*, "response")= int 1
  .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. .. ..- attr(*, "predvars")= language list(YIELD, center(FERTILIZER, center = 137.5, verbose = FALSE))
  .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
  .. .. .. ..- attr(*, "names")= chr [1:2] "YIELD" "center(FERTILIZER)"
 - attr(*, "class")= chr "lm"
## Return the data used to fit the model 
fert_lm1$model
   YIELD center(FERTILIZER)
1     84             -112.5
2     80              -87.5
3     90              -62.5
4    154              -37.5
5    148              -12.5
6    169               12.5
7    206               37.5
8    244               62.5
9    212               87.5
10   248              112.5
## Return the estimated coefficients
fert_lm1$coefficients
       (Intercept) center(FERTILIZER) 
       163.5000000          0.8113939 

Rather than directly access these attributes (which may have different names in different modelling functions), there are extraction functions that reach into the model output and extract the information in a convenient way.

## Extract the estimated coefficients
coef(fert_lm1)
       (Intercept) center(FERTILIZER) 
       163.5000000          0.8113939 
## Extract the fitted (predicted) values
fitted(fert_lm1)
        1         2         3         4         5         6         7         8 
 72.21818  92.50303 112.78788 133.07273 153.35758 173.64242 193.92727 214.21212 
        9        10 
234.49697 254.78182 
## Extract the residules
resid(fert_lm1)
         1          2          3          4          5          6          7 
 11.781818 -12.503030 -22.787879  20.927273  -5.357576  -4.642424  12.072727 
         8          9         10 
 29.787879 -22.496970  -6.781818 
args(residuals.lm)
function (object, type = c("working", "response", "deviance", 
    "pearson", "partial"), ...) 
NULL
fert_glmmTMB <- glmmTMB::glmmTMB(YIELD ~ 1 + FERTILIZER, data = fert)
fert_glmmTMB <- glmmTMB::glmmTMB(YIELD ~ FERTILIZER, data = fert)

The glmmTMB() function returns a list comprising the estimated coefficients, the residuals, the data used to fit the model and various other attributes associated with the model fitting procedure.

## Get the name of the attributes within the glmmTMB model
attributes(fert_glmmTMB)
$names
[1] "obj"       "fit"       "sdr"       "call"      "frame"     "modelInfo"
[7] "fitted"   

$class
[1] "glmmTMB"
## Explore all the conttents of the fitted model
str(fert_glmmTMB)
List of 7
 $ obj      :List of 10
  ..$ par     : Named num [1:3] 1 1 0
  .. ..- attr(*, "names")= chr [1:3] "beta" "beta" "betad"
  ..$ fn      :function (x = last.par[lfixed()], ...)  
  ..$ gr      :function (x = last.par[lfixed()], ...)  
  ..$ he      :function (x = last.par[lfixed()], atomic = usingAtomics())  
  ..$ hessian : logi FALSE
  ..$ method  : chr "BFGS"
  ..$ retape  :function (set.defaults = TRUE)  
  ..$ env     :<environment: 0x563bf2ee89c0> 
  ..$ report  :function (par = last.par)  
  ..$ simulate:function (par = last.par, complete = FALSE)  
 $ fit      :List of 7
  ..$ par        : Named num [1:3] 51.933 0.811 5.666
  .. ..- attr(*, "names")= chr [1:3] "beta" "beta" "betad"
  ..$ objective  : num 42.5
  ..$ convergence: int 0
  ..$ iterations : int 33
  ..$ evaluations: Named int [1:2] 37 34
  .. ..- attr(*, "names")= chr [1:2] "function" "gradient"
  ..$ message    : chr "relative convergence (4)"
  ..$ parfull    : Named num [1:3] 51.933 0.811 5.666
  .. ..- attr(*, "names")= chr [1:3] "beta" "beta" "betad"
 $ sdr      :List of 8
  ..$ value         : num(0) 
  ..$ sd            : num(0) 
  ..$ cov           : logi[0 , 0 ] 
  ..$ par.fixed     : Named num [1:3] 51.933 0.811 5.666
  .. ..- attr(*, "names")= chr [1:3] "beta" "beta" "betad"
  ..$ cov.fixed     : num [1:3, 1:3] 1.35e+02 -7.70e-01 -3.53e-06 -7.70e-01 5.60e-03 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:3] "beta" "beta" "betad"
  .. .. ..$ : chr [1:3] "beta" "beta" "betad"
  ..$ pdHess        : logi TRUE
  ..$ gradient.fixed: num [1, 1:3] -1.88e-07 -9.95e-06 -3.76e-07
  ..$ env           :<environment: 0x563bf2849890> 
  ..- attr(*, "class")= chr "sdreport"
 $ call     : language glmmTMB::glmmTMB(formula = YIELD ~ FERTILIZER, data = fert, ziformula = ~0,      dispformula = ~1)
 $ frame    :'data.frame':  10 obs. of  2 variables:
  ..$ YIELD     : num [1:10] 84 80 90 154 148 169 206 244 212 248
  ..$ FERTILIZER: num [1:10] 25 50 75 100 125 150 175 200 225 250
  ..- attr(*, "terms")=Classes 'terms', 'formula'  language YIELD ~ FERTILIZER + 0 + 1
  .. .. ..- attr(*, "variables")= language list(YIELD, FERTILIZER)
  .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:2] "YIELD" "FERTILIZER"
  .. .. .. .. ..$ : chr "FERTILIZER"
  .. .. ..- attr(*, "term.labels")= chr "FERTILIZER"
  .. .. ..- attr(*, "order")= int 1
  .. .. ..- attr(*, "intercept")= int 1
  .. .. ..- attr(*, "response")= int 1
  .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. .. ..- attr(*, "predvars")= language list(YIELD, FERTILIZER)
  .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
  .. .. .. ..- attr(*, "names")= chr [1:2] "YIELD" "FERTILIZER"
 $ modelInfo:List of 13
  ..$ nobs     : int 10
  ..$ respCol  : Named int 1
  .. ..- attr(*, "names")= chr "YIELD"
  ..$ grpVar   : NULL
  ..$ family   :List of 12
  .. ..$ family    : chr "gaussian"
  .. ..$ link      : chr "identity"
  .. ..$ linkfun   :function (mu)  
  .. ..$ linkinv   :function (eta)  
  .. ..$ variance  :function (mu)  
  .. ..$ dev.resids:function (y, mu, wt)  
  .. ..$ aic       :function (y, n, mu, wt, dev)  
  .. ..$ mu.eta    :function (eta)  
  .. ..$ initialize:  expression({  n <- rep.int(1, nobs)  if (is.null(etastart) && is.null(start) && is.null(mustart) && ((family$link| __truncated__
  .. ..$ validmu   :function (mu)  
  .. ..$ valideta  :function (eta)  
  .. ..$ dispersion: num NA
  .. ..- attr(*, "class")= chr "family"
  ..$ contrasts: NULL
  ..$ reTrms   :List of 2
  .. ..$ cond:List of 1
  .. .. ..$ terms:List of 1
  .. .. .. ..$ fixed:Classes 'terms', 'formula'  language YIELD ~ FERTILIZER
  .. .. .. .. .. ..- attr(*, "variables")= language list(YIELD, FERTILIZER)
  .. .. .. .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. .. .. ..$ : chr [1:2] "YIELD" "FERTILIZER"
  .. .. .. .. .. .. .. ..$ : chr "FERTILIZER"
  .. .. .. .. .. ..- attr(*, "term.labels")= chr "FERTILIZER"
  .. .. .. .. .. ..- attr(*, "order")= int 1
  .. .. .. .. .. ..- attr(*, "intercept")= int 1
  .. .. .. .. .. ..- attr(*, "response")= int 1
  .. .. .. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. .. .. .. .. ..- attr(*, "predvars")= language list(YIELD, FERTILIZER)
  .. .. .. .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
  .. .. .. .. .. .. ..- attr(*, "names")= chr [1:2] "YIELD" "FERTILIZER"
  .. ..$ zi  :List of 1
  .. .. ..$ terms: NULL
  ..$ terms    :List of 3
  .. ..$ cond:List of 1
  .. .. ..$ fixed:Classes 'terms', 'formula'  language YIELD ~ FERTILIZER
  .. .. .. .. ..- attr(*, "variables")= language list(YIELD, FERTILIZER)
  .. .. .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. .. ..$ : chr [1:2] "YIELD" "FERTILIZER"
  .. .. .. .. .. .. ..$ : chr "FERTILIZER"
  .. .. .. .. ..- attr(*, "term.labels")= chr "FERTILIZER"
  .. .. .. .. ..- attr(*, "order")= int 1
  .. .. .. .. ..- attr(*, "intercept")= int 1
  .. .. .. .. ..- attr(*, "response")= int 1
  .. .. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. .. .. .. ..- attr(*, "predvars")= language list(YIELD, FERTILIZER)
  .. .. .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
  .. .. .. .. .. ..- attr(*, "names")= chr [1:2] "YIELD" "FERTILIZER"
  .. ..$ zi  : NULL
  .. ..$ disp:List of 1
  .. .. ..$ fixed:Classes 'terms', 'formula'  language ~1
  .. .. .. .. ..- attr(*, "variables")= language list()
  .. .. .. .. ..- attr(*, "factors")= int(0) 
  .. .. .. .. ..- attr(*, "term.labels")= chr(0) 
  .. .. .. .. ..- attr(*, "order")= int(0) 
  .. .. .. .. ..- attr(*, "intercept")= int 1
  .. .. .. .. ..- attr(*, "response")= int 0
  .. .. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. .. .. .. ..- attr(*, "predvars")= language list()
  .. .. .. .. ..- attr(*, "dataClasses")= Named chr(0) 
  .. .. .. .. .. ..- attr(*, "names")= chr(0) 
  ..$ reStruc  :List of 2
  .. ..$ condReStruc: list()
  .. ..$ ziReStruc  : list()
  ..$ allForm  :List of 4
  .. ..$ combForm   :Class 'formula'  language YIELD ~ FERTILIZER + 0 + 1
  .. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. ..$ formula    :Class 'formula'  language YIELD ~ FERTILIZER
  .. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. ..$ ziformula  :Class 'formula'  language ~0
  .. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. ..$ dispformula:Class 'formula'  language ~1
  .. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  ..$ REML     : logi FALSE
  ..$ map      : NULL
  ..$ sparseX  : Named logi [1:3] FALSE FALSE FALSE
  .. ..- attr(*, "names")= chr [1:3] "cond" "zi" "disp"
  ..$ parallel : int 1
 $ fitted   : NULL
 - attr(*, "class")= chr "glmmTMB"
## Return the data used to fit the model 
fert_glmmTMB$frame
   YIELD FERTILIZER
1     84         25
2     80         50
3     90         75
4    154        100
5    148        125
6    169        150
7    206        175
8    244        200
9    212        225
10   248        250
## Return the estimated coefficients
fert_glmmTMB$fit
$par
     beta      beta     betad 
51.933316  0.811394  5.665667 

$objective
[1] 42.51772

$convergence
[1] 0

$iterations
[1] 33

$evaluations
function gradient 
      37       34 

$message
[1] "relative convergence (4)"

$parfull
     beta      beta     betad 
51.933316  0.811394  5.665667 

Rather than directly access these attributes (which may have different names in different modelling functions), there are extraction functions that reach into the model output and extract the information in a convenient way.

## Extract the estimated coefficients
fixef(fert_glmmTMB)

Conditional model:
(Intercept)   FERTILIZER  
    51.9333       0.8114  
## Extract the fitted (predicted) values
fitted(fert_glmmTMB)
 [1]  72.21817  92.50302 112.78787 133.07272 153.35757 173.64242 193.92727
 [8] 214.21212 234.49697 254.78182
## Extract the residules
resid(fert_glmmTMB)
         1          2          3          4          5          6          7 
 11.781834 -12.503017 -22.787868  20.927281  -5.357569  -4.642420  12.072729 
         8          9         10 
 29.787879 -22.496972  -6.781823 
args(glmmTMB:::residuals.glmmTMB)
function (object, type = c("response", "pearson", "working"), 
    ...) 
NULL
fert_glmmTMB1 <- glmmTMB::glmmTMB(YIELD ~ 1 + center(FERTILIZER), data = fert)
fert_glmmTMB1 <- glmmTMB::glmmTMB(YIELD ~ center(FERTILIZER), data = fert)

The glmmTMB() function returns a list comprising the estimated coefficients, the residuals, the data used to fit the model and various other attributes associated with the model fitting procedure.

## Get the name of the attributes within the glmmTMB model
attributes(fert_glmmTMB1)
$names
[1] "obj"       "fit"       "sdr"       "call"      "frame"     "modelInfo"
[7] "fitted"   

$class
[1] "glmmTMB"
## Explore all the conttents of the fitted model
str(fert_glmmTMB1)
List of 7
 $ obj      :List of 10
  ..$ par     : Named num [1:3] 1 1 0
  .. ..- attr(*, "names")= chr [1:3] "beta" "beta" "betad"
  ..$ fn      :function (x = last.par[lfixed()], ...)  
  ..$ gr      :function (x = last.par[lfixed()], ...)  
  ..$ he      :function (x = last.par[lfixed()], atomic = usingAtomics())  
  ..$ hessian : logi FALSE
  ..$ method  : chr "BFGS"
  ..$ retape  :function (set.defaults = TRUE)  
  ..$ env     :<environment: 0x555d629c1218> 
  ..$ report  :function (par = last.par)  
  ..$ simulate:function (par = last.par, complete = FALSE)  
 $ fit      :List of 7
  ..$ par        : Named num [1:3] 163.5 0.811 5.666
  .. ..- attr(*, "names")= chr [1:3] "beta" "beta" "betad"
  ..$ objective  : num 42.5
  ..$ convergence: int 0
  ..$ iterations : int 46
  ..$ evaluations: Named int [1:2] 54 47
  .. ..- attr(*, "names")= chr [1:2] "function" "gradient"
  ..$ message    : chr "relative convergence (4)"
  ..$ parfull    : Named num [1:3] 163.5 0.811 5.666
  .. ..- attr(*, "names")= chr [1:3] "beta" "beta" "betad"
 $ sdr      :List of 8
  ..$ value         : num(0) 
  ..$ sd            : num(0) 
  ..$ cov           : logi[0 , 0 ] 
  ..$ par.fixed     : Named num [1:3] 163.5 0.811 5.666
  .. ..- attr(*, "names")= chr [1:3] "beta" "beta" "betad"
  ..$ cov.fixed     : num [1:3, 1:3] 2.89e+01 1.46e-13 1.67e-07 1.46e-13 5.60e-03 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:3] "beta" "beta" "betad"
  .. .. ..$ : chr [1:3] "beta" "beta" "betad"
  ..$ pdHess        : logi TRUE
  ..$ gradient.fixed: num [1, 1:3] 2.88e-08 3.56e-06 -3.31e-07
  ..$ env           :<environment: 0x555d639e9530> 
  ..- attr(*, "class")= chr "sdreport"
 $ call     : language glmmTMB::glmmTMB(formula = YIELD ~ center(FERTILIZER), data = fert, ziformula = ~0,      dispformula = ~1)
 $ frame    :'data.frame':  10 obs. of  2 variables:
  ..$ YIELD             : num [1:10] 84 80 90 154 148 169 206 244 212 248
  ..$ center(FERTILIZER): 'dw_transformer' num [1:10] -112.5 -87.5 -62.5 -37.5 -12.5 ...
  .. ..- attr(*, "center")= num 138
  .. ..- attr(*, "scale")= num 1
  .. ..- attr(*, "robust")= logi FALSE
  ..- attr(*, "terms")=Classes 'terms', 'formula'  language YIELD ~ center(FERTILIZER) + 0 + 1
  .. .. ..- attr(*, "variables")= language list(YIELD, center(FERTILIZER))
  .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. ..$ : chr [1:2] "YIELD" "center(FERTILIZER)"
  .. .. .. .. ..$ : chr "center(FERTILIZER)"
  .. .. ..- attr(*, "term.labels")= chr "center(FERTILIZER)"
  .. .. ..- attr(*, "order")= int 1
  .. .. ..- attr(*, "intercept")= int 1
  .. .. ..- attr(*, "response")= int 1
  .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. .. ..- attr(*, "predvars")= language list(YIELD, center(FERTILIZER, center = 137.5, verbose = FALSE))
  .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
  .. .. .. ..- attr(*, "names")= chr [1:2] "YIELD" "center(FERTILIZER)"
 $ modelInfo:List of 13
  ..$ nobs     : int 10
  ..$ respCol  : Named int 1
  .. ..- attr(*, "names")= chr "YIELD"
  ..$ grpVar   : NULL
  ..$ family   :List of 12
  .. ..$ family    : chr "gaussian"
  .. ..$ link      : chr "identity"
  .. ..$ linkfun   :function (mu)  
  .. ..$ linkinv   :function (eta)  
  .. ..$ variance  :function (mu)  
  .. ..$ dev.resids:function (y, mu, wt)  
  .. ..$ aic       :function (y, n, mu, wt, dev)  
  .. ..$ mu.eta    :function (eta)  
  .. ..$ initialize:  expression({  n <- rep.int(1, nobs)  if (is.null(etastart) && is.null(start) && is.null(mustart) && ((family$link| __truncated__
  .. ..$ validmu   :function (mu)  
  .. ..$ valideta  :function (eta)  
  .. ..$ dispersion: num NA
  .. ..- attr(*, "class")= chr "family"
  ..$ contrasts: NULL
  ..$ reTrms   :List of 2
  .. ..$ cond:List of 1
  .. .. ..$ terms:List of 1
  .. .. .. ..$ fixed:Classes 'terms', 'formula'  language YIELD ~ center(FERTILIZER)
  .. .. .. .. .. ..- attr(*, "variables")= language list(YIELD, center(FERTILIZER))
  .. .. .. .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. .. .. ..$ : chr [1:2] "YIELD" "center(FERTILIZER)"
  .. .. .. .. .. .. .. ..$ : chr "center(FERTILIZER)"
  .. .. .. .. .. ..- attr(*, "term.labels")= chr "center(FERTILIZER)"
  .. .. .. .. .. ..- attr(*, "order")= int 1
  .. .. .. .. .. ..- attr(*, "intercept")= int 1
  .. .. .. .. .. ..- attr(*, "response")= int 1
  .. .. .. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. .. .. .. .. ..- attr(*, "predvars")= language list(YIELD, center(FERTILIZER, center = 137.5, verbose = FALSE))
  .. .. .. .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
  .. .. .. .. .. .. ..- attr(*, "names")= chr [1:2] "YIELD" "center(FERTILIZER)"
  .. ..$ zi  :List of 1
  .. .. ..$ terms: NULL
  ..$ terms    :List of 3
  .. ..$ cond:List of 1
  .. .. ..$ fixed:Classes 'terms', 'formula'  language YIELD ~ center(FERTILIZER)
  .. .. .. .. ..- attr(*, "variables")= language list(YIELD, center(FERTILIZER))
  .. .. .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
  .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. .. .. .. .. ..$ : chr [1:2] "YIELD" "center(FERTILIZER)"
  .. .. .. .. .. .. ..$ : chr "center(FERTILIZER)"
  .. .. .. .. ..- attr(*, "term.labels")= chr "center(FERTILIZER)"
  .. .. .. .. ..- attr(*, "order")= int 1
  .. .. .. .. ..- attr(*, "intercept")= int 1
  .. .. .. .. ..- attr(*, "response")= int 1
  .. .. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. .. .. .. ..- attr(*, "predvars")= language list(YIELD, center(FERTILIZER, center = 137.5, verbose = FALSE))
  .. .. .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
  .. .. .. .. .. ..- attr(*, "names")= chr [1:2] "YIELD" "center(FERTILIZER)"
  .. ..$ zi  : NULL
  .. ..$ disp:List of 1
  .. .. ..$ fixed:Classes 'terms', 'formula'  language ~1
  .. .. .. .. ..- attr(*, "variables")= language list()
  .. .. .. .. ..- attr(*, "factors")= int(0) 
  .. .. .. .. ..- attr(*, "term.labels")= chr(0) 
  .. .. .. .. ..- attr(*, "order")= int(0) 
  .. .. .. .. ..- attr(*, "intercept")= int 1
  .. .. .. .. ..- attr(*, "response")= int 0
  .. .. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. .. .. .. ..- attr(*, "predvars")= language list()
  .. .. .. .. ..- attr(*, "dataClasses")= Named chr(0) 
  .. .. .. .. .. ..- attr(*, "names")= chr(0) 
  ..$ reStruc  :List of 2
  .. ..$ condReStruc: list()
  .. ..$ ziReStruc  : list()
  ..$ allForm  :List of 4
  .. ..$ combForm   :Class 'formula'  language YIELD ~ center(FERTILIZER) + 0 + 1
  .. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. ..$ formula    :Class 'formula'  language YIELD ~ center(FERTILIZER)
  .. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. ..$ ziformula  :Class 'formula'  language ~0
  .. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  .. ..$ dispformula:Class 'formula'  language ~1
  .. .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv> 
  ..$ REML     : logi FALSE
  ..$ map      : NULL
  ..$ sparseX  : Named logi [1:3] FALSE FALSE FALSE
  .. ..- attr(*, "names")= chr [1:3] "cond" "zi" "disp"
  ..$ parallel : int 1
 $ fitted   : NULL
 - attr(*, "class")= chr "glmmTMB"
## Return the data used to fit the model 
fert_glmmTMB1$frame
   YIELD center(FERTILIZER)
1     84             -112.5
2     80              -87.5
3     90              -62.5
4    154              -37.5
5    148              -12.5
6    169               12.5
7    206               37.5
8    244               62.5
9    212               87.5
10   248              112.5
## Return the estimated coefficients
fert_glmmTMB1$fit
$par
      beta       beta      betad 
163.500001   0.811394   5.665667 

$objective
[1] 42.51772

$convergence
[1] 0

$iterations
[1] 46

$evaluations
function gradient 
      54       47 

$message
[1] "relative convergence (4)"

$parfull
      beta       beta      betad 
163.500001   0.811394   5.665667 

Rather than directly access these attributes (which may have different names in different modelling functions), there are extraction functions that reach into the model output and extract the information in a convenient way.

## Extract the estimated coefficients
fixef(fert_glmmTMB1)

Conditional model:
       (Intercept)  center(FERTILIZER)  
          163.5000              0.8114  
## Extract the fitted (predicted) values
fitted(fert_glmmTMB1)
 [1]  72.21818  92.50303 112.78788 133.07273 153.35758 173.64243 193.92727
 [8] 214.21212 234.49697 254.78182
## Extract the residules
resid(fert_glmmTMB1)
         1          2          3          4          5          6          7 
 11.781820 -12.503029 -22.787878  20.927273  -5.357576  -4.642425  12.072726 
         8          9         10 
 29.787877 -22.496972  -6.781821 
args(glmmTMB:::residuals.glmmTMB)
function (object, type = c("response", "pearson", "working"), 
    ...) 
NULL

Lets performing Ordinary Least Squares (OLS) regression by first principles.

To do so, we expresses the response as a vector (\(Y\)), the deterministic (linear predictor) as a matrix (\(\boldsymbol{X}\)) comprising of a column of 1’s (for multiplying the intercept parameter) and the predictor variable (for multiplying the slope parameter). We then also have a vector of beta parameters (\(\beta_0\) and \(\beta\) representing the intercept and slope respectively).

\[ y_i = \beta_0\times 1 + \beta_1\times x_i \]

\[ \underbrace{\begin{bmatrix} y_1\\ y_2\\ y_3\\ ...\\ y_i \end{bmatrix}}_{Y} = \underbrace{\begin{bmatrix} 1&x_1\\ 1&x_2\\ 1&x_3\\ ...\\ 1&x_i \end{bmatrix}}_{\boldsymbol{X}} \underbrace{\vphantom{\begin{bmatrix} y_1\\ y_2\\ y_3\\ ...\\ y_i \end{bmatrix}}\begin{bmatrix} \beta_0&\beta \end{bmatrix}}_{\boldsymbol{\beta}} \]

In OLS we estimate the parameters (\(\beta_0\) and \(\beta\)) by minimising the sum of the squared residuals (\(\boldsymbol{\varepsilon}\)).

\[ \begin{align} Y =& \boldsymbol{X}\boldsymbol{\beta} + \boldsymbol{\varepsilon}\\ \boldsymbol{\varepsilon} =& Y - \boldsymbol{X}\boldsymbol{\beta} \end{align} \]

We can then minimise the sum of squares residuals. This is essentially \(\boldsymbol{\varepsilon}^T\boldsymbol{\varepsilon}\), where \(\boldsymbol{\varepsilon}\) is \(\boldsymbol{\varepsilon}\) transposed.

\[ \small \begin{bmatrix} \varepsilon_1&\varepsilon_2&\varepsilon_3&...&\varepsilon_n \end{bmatrix} \begin{bmatrix} \varepsilon_1\\ \varepsilon_2\\ \varepsilon_3\\ ...\\ \varepsilon_n\\ \end{bmatrix} = \begin{bmatrix} \varepsilon_1\times\varepsilon_1 + \varepsilon_2\times\varepsilon_2 + \varepsilon_3\times\varepsilon_3 + ... + \varepsilon_n\times\varepsilon_n \end{bmatrix} \]

\[ \begin{align} \boldsymbol{\varepsilon}^T\boldsymbol{\varepsilon} =& (Y - \boldsymbol{X}\boldsymbol{\beta})^T(Y - \boldsymbol{X}\boldsymbol{\beta})\\ \boldsymbol{\beta} =& (\boldsymbol{X}^T\boldsymbol{X})^{-1} \boldsymbol{X}^T Y \end{align} \]

## Generate the model matrix
X <- model.matrix(~FERTILIZER, data = fert)
## Solve for beta
beta <- solve(t(X) %*% X) %*% t(X) %*% fert$YIELD
beta
                  [,1]
(Intercept) 51.9333333
FERTILIZER   0.8113939

The above simple linear regression model can be fit using the lm() function. This function, uses Ordinary Least Squares (OLS). The full model formula would include a 1 (for the intercept) and the predictor values (for the slope). In R, models generally always include an intercept by default and thus we can abbreviate the formula a little by omitting the 1.

6 Model validation

7 Model outputs

Prior to examining the estimated coefficients and any associated hypothesis tests, it is a good idea to explore the fitted trends. Not only does this give you a sense for the overall outcomes (and help you interpret the model summaries etc), it also provides an opportunity to reflect on whether the model has yielded sensible patterns.

There are numerous ways to explore the fitted trends and these are largely based on two types of effects:

  • conditional effects: predictions that are conditioned on certain levels (typically the reference level) of factors. For example, the trend/effects of one predictor at the first level (or first combination) of other categorical predictor(s).
  • marginal effects: predictions that are marginalised (= averaged) over all levels of the factors. For example, the trend/effects of one predictor averaged across all levels (or combinations) of other predictors.

There are also numerous routines and packages to support these fitted trends (partial plots).

Package Function Type Notes
sjPlot plot_model(type = 'eff') Marginal means
effects allEffects() Marginal means
ggeffects ggeffects() Marginal means calls effects()
ggeffects ggpredict() Conditional means calls predict()
ggeffects ggemmeans() Marginal means calls emmeans()

Note, in the absence of categorical predictors, they will all yield the same…

8 Model investigation / hypothesis testing

9 Predictions

For simple models prediction is essentially taking the model formula complete with parameter (coefficient) estimates and solving for new values of the predictor. To explore this, we will use the fitted model to predict Yield for a Fertiliser concentration of 110.

10 Additional analyses

11 Summary figures

Although there are numerous easy to use routines that will generate partial plots (see above), it is also useful to be able to produce the data behind the figures yourself. That way, you can have more control over the type and style of the figures.

Producing a summary figure is essentially plotting the results of predictions. In the case of a trend, we want a series of predictions associated with a sequence of predictor values. Hence, we establish a prediction grid that contains the sequence of values for the predictor we would like to display.

There are a number of functions we can use to generate different sequences of prediction values. For example, we may want a sequence of values from the lowest to the highest observed predictor value. Alternatively, we may just want to explore predictions at the first, third and fifth quantiles. The following table indicates some of the functions that are helpful for these purposes.

Function Values returned
modelr::seq_range() equal increment values from smallest to largest
Hmisc::smean_sdl() mean as well as mean plus/minus one standard deviation
Hmisc::smedian_hilow() median as well as min and max
Hmisc::smean.cl.normal() mean as well as lower and upper 95% confidence interval

12 Reporting

report::report(fert_lm)
We fitted a linear model (estimated using OLS) to predict YIELD with FERTILIZER
(formula: YIELD ~ FERTILIZER). The model explains a statistically significant
and substantial proportion of variance (R2 = 0.92, F(1, 8) = 94.04, p < .001,
adj. R2 = 0.91). The model's intercept, corresponding to FERTILIZER = 0, is at
51.93 (95% CI [22.00, 81.86], t(8) = 4.00, p = 0.004). Within this model:

  - The effect of FERTILIZER is statistically significant and positive (beta =
0.81, 95% CI [0.62, 1.00], t(8) = 9.70, p < .001; Std. beta = 0.96, 95% CI
[0.73, 1.19])

Standardized parameters were obtained by fitting the model on a standardized
version of the dataset. 95% Confidence Intervals (CIs) and p-values were
computed using a Wald t-distribution approximation.
report::report_text(fert_lm)
We fitted a linear model (estimated using OLS) to predict YIELD with FERTILIZER
(formula: YIELD ~ FERTILIZER). The model explains a statistically significant
and substantial proportion of variance (R2 = 0.92, F(1, 8) = 94.04, p < .001,
adj. R2 = 0.91). The model's intercept, corresponding to FERTILIZER = 0, is at
51.93 (95% CI [22.00, 81.86], t(8) = 4.00, p = 0.004). Within this model:

  - The effect of FERTILIZER is statistically significant and positive (beta =
0.81, 95% CI [0.62, 1.00], t(8) = 9.70, p < .001; Std. beta = 0.96, 95% CI
[0.73, 1.19])

Standardized parameters were obtained by fitting the model on a standardized
version of the dataset. 95% Confidence Intervals (CIs) and p-values were
computed using a Wald t-distribution approximation.
report::report_model(fert_lm)
linear model (estimated using OLS) to predict YIELD with FERTILIZER (formula: YIELD ~ FERTILIZER)
report::report_parameters(fert_lm)
  - The intercept is statistically significant and positive (beta = 51.93, 95% CI [22.00, 81.86], t(8) = 4.00, p = 0.004; Std. beta = 0.00, 95% CI [-0.22, 0.22])
  - The effect of FERTILIZER is statistically significant and positive (beta = 0.81, 95% CI [0.62, 1.00], t(8) = 9.70, p < .001; Std. beta = 0.96, 95% CI [0.73, 1.19])
report::report_statistics(fert_lm)
beta = 51.93, 95% CI [22.00, 81.86], t(8) = 4.00, p = 0.004; Std. beta = 0.00, 95% CI [-0.22, 0.22]
beta = 0.81, 95% CI [0.62, 1.00], t(8) = 9.70, p < .001; Std. beta = 0.96, 95% CI [0.73, 1.19]
report::report_info(fert_lm)
Standardized parameters were obtained by fitting the model on a standardized version of the dataset. 95% Confidence Intervals (CIs) and p-values were computed using a Wald t-distribution approximation.
report::report_performance(fert_lm)
The model explains a statistically significant and substantial proportion of
variance (R2 = 0.92, F(1, 8) = 94.04, p < .001, adj. R2 = 0.91)
report::report_system()
Analyses were conducted using the R Statistical language (version 4.3.1; R Core
Team, 2023) on ArcoLinux
report::report(sessionInfo())
Analyses were conducted using the R Statistical language (version 4.3.1; R Core
Team, 2023) on ArcoLinux, using the packages effectsize (version 0.8.3;
Ben-Shachar MS et al., 2020), glmmTMB (version 1.1.7; Brooks ME et al., 2017),
effects (version 4.2.2; Fox J, Weisberg S, 2019), car (version 3.1.2; Fox J,
Weisberg S, 2019), carData (version 3.0.5; Fox J et al., 2022), lubridate
(version 1.9.2; Grolemund G, Wickham H, 2011), DHARMa (version 0.4.6; Hartig F,
2022), lindia (version 0.9; Lee Y, Ventura S, 2017), emmeans (version 1.8.7;
Lenth R, 2023), ggeffects (version 1.2.3; Lüdecke D, 2018), sjPlot (version
2.8.14; Lüdecke D, 2023), parameters (version 0.21.1; Lüdecke D et al., 2020),
performance (version 0.10.4; Lüdecke D et al., 2021), easystats (version 0.6.0;
Lüdecke D et al., 2022), see (version 0.7.4; Lüdecke D et al., 2021), insight
(version 0.19.3; Lüdecke D et al., 2019), bayestestR (version 0.13.0; Makowski
D et al., 2019), modelbased (version 0.8.6; Makowski D et al., 2020),
correlation (version 0.8.3; Makowski D et al., 2020), report (version 0.5.7;
Makowski D et al., 2023), tibble (version 3.2.1; Müller K, Wickham H, 2023),
datawizard (version 0.8.0; Patil I et al., 2022), broom (version 1.0.5;
Robinson D et al., 2023), ggfortify (version 0.4.16; Tang Y et al., 2016),
ggplot2 (version 3.4.2; Wickham H, 2016), stringr (version 1.5.0; Wickham H,
2022), forcats (version 1.0.0; Wickham H, 2023), modelr (version 0.1.11;
Wickham H, 2023), tidyverse (version 2.0.0; Wickham H et al., 2019), dplyr
(version 1.1.2; Wickham H et al., 2023), purrr (version 1.0.1; Wickham H, Henry
L, 2023), readr (version 2.1.4; Wickham H et al., 2023), tidyr (version 1.3.0;
Wickham H et al., 2023) and knitr (version 1.43; Xie Y, 2023).

References
----------
  - Ben-Shachar MS, Lüdecke D, Makowski D (2020). "effectsize: Estimation of
Effect Size Indices and Standardized Parameters." _Journal of Open Source
Software_, *5*(56), 2815. doi:10.21105/joss.02815
<https://doi.org/10.21105/joss.02815>, <https://doi.org/10.21105/joss.02815>.
  - Brooks ME, Kristensen K, van Benthem KJ, Magnusson A, Berg CW, Nielsen A,
Skaug HJ, Maechler M, Bolker BM (2017). "glmmTMB Balances Speed and Flexibility
Among Packages for Zero-inflated Generalized Linear Mixed Modeling." _The R
Journal_, *9*(2), 378-400. doi:10.32614/RJ-2017-066
<https://doi.org/10.32614/RJ-2017-066>.
  - Fox J, Weisberg S (2019). _An R Companion to Applied Regression_, 3rd
edition. Sage, Thousand Oaks CA.
<https://socialsciences.mcmaster.ca/jfox/Books/Companion/index.html>. Fox J,
Weisberg S (2018). "Visualizing Fit and Lack of Fit in Complex Regression
Models with Predictor Effect Plots and Partial Residuals." _Journal of
Statistical Software_, *87*(9), 1-27. doi:10.18637/jss.v087.i09
<https://doi.org/10.18637/jss.v087.i09>. Fox J (2003). "Effect Displays in R
for Generalised Linear Models." _Journal of Statistical Software_, *8*(15),
1-27. doi:10.18637/jss.v008.i15 <https://doi.org/10.18637/jss.v008.i15>. Fox J,
Hong J (2009). "Effect Displays in R for Multinomial and Proportional-Odds
Logit Models: Extensions to the effects Package." _Journal of Statistical
Software_, *32*(1), 1-24. doi:10.18637/jss.v032.i01
<https://doi.org/10.18637/jss.v032.i01>.
  - Fox J, Weisberg S (2019). _An R Companion to Applied Regression_, Third
edition. Sage, Thousand Oaks CA.
<https://socialsciences.mcmaster.ca/jfox/Books/Companion/>.
  - Fox J, Weisberg S, Price B (2022). _carData: Companion to Applied Regression
Data Sets_. R package version 3.0-5,
<https://CRAN.R-project.org/package=carData>.
  - Grolemund G, Wickham H (2011). "Dates and Times Made Easy with lubridate."
_Journal of Statistical Software_, *40*(3), 1-25.
<https://www.jstatsoft.org/v40/i03/>.
  - Hartig F (2022). _DHARMa: Residual Diagnostics for Hierarchical (Multi-Level
/ Mixed) Regression Models_. R package version 0.4.6,
<https://CRAN.R-project.org/package=DHARMa>.
  - Lee Y, Ventura S (2017). _lindia: Automated Linear Regression Diagnostic_. R
package version 0.9, <https://CRAN.R-project.org/package=lindia>.
  - Lenth R (2023). _emmeans: Estimated Marginal Means, aka Least-Squares Means_.
R package version 1.8.7, <https://CRAN.R-project.org/package=emmeans>.
  - Lüdecke D (2018). "ggeffects: Tidy Data Frames of Marginal Effects from
Regression Models." _Journal of Open Source Software_, *3*(26), 772.
doi:10.21105/joss.00772 <https://doi.org/10.21105/joss.00772>.
  - Lüdecke D (2023). _sjPlot: Data Visualization for Statistics in Social
Science_. R package version 2.8.14,
<https://CRAN.R-project.org/package=sjPlot>.
  - Lüdecke D, Ben-Shachar M, Patil I, Makowski D (2020). "Extracting, Computing
and Exploring the Parameters of Statistical Models using R." _Journal of Open
Source Software_, *5*(53), 2445. doi:10.21105/joss.02445
<https://doi.org/10.21105/joss.02445>.
  - Lüdecke D, Ben-Shachar M, Patil I, Waggoner P, Makowski D (2021).
"performance: An R Package for Assessment, Comparison and Testing of
Statistical Models." _Journal of Open Source Software_, *6*(60), 3139.
doi:10.21105/joss.03139 <https://doi.org/10.21105/joss.03139>.
  - Lüdecke D, Ben-Shachar M, Patil I, Wiernik B, Makowski D (2022). "easystats:
Framework for Easy Statistical Modeling, Visualization, and Reporting." _CRAN_.
R package, <https://easystats.github.io/easystats/>.
  - Lüdecke D, Patil I, Ben-Shachar M, Wiernik B, Waggoner P, Makowski D (2021).
"see: An R Package for Visualizing Statistical Models." _Journal of Open Source
Software_, *6*(64), 3393. doi:10.21105/joss.03393
<https://doi.org/10.21105/joss.03393>.
  - Lüdecke D, Waggoner P, Makowski D (2019). "insight: A Unified Interface to
Access Information from Model Objects in R." _Journal of Open Source Software_,
*4*(38), 1412. doi:10.21105/joss.01412 <https://doi.org/10.21105/joss.01412>.
  - Makowski D, Ben-Shachar M, Lüdecke D (2019). "bayestestR: Describing Effects
and their Uncertainty, Existence and Significance within the Bayesian
Framework." _Journal of Open Source Software_, *4*(40), 1541.
doi:10.21105/joss.01541 <https://doi.org/10.21105/joss.01541>,
<https://joss.theoj.org/papers/10.21105/joss.01541>.
  - Makowski D, Ben-Shachar M, Patil I, Lüdecke D (2020). "Estimation of
Model-Based Predictions, Contrasts and Means." _CRAN_.
<https://github.com/easystats/modelbased>.
  - Makowski D, Ben-Shachar M, Patil I, Lüdecke D (2020). "Methods and Algorithms
for Correlation Analysis in R." _Journal of Open Source Software_, *5*(51),
2306. doi:10.21105/joss.02306 <https://doi.org/10.21105/joss.02306>,
<https://joss.theoj.org/papers/10.21105/joss.02306>.
  - Makowski D, Lüdecke D, Patil I, Thériault R, Ben-Shachar M, Wiernik B (2023).
"Automated Results Reporting as a Practical Tool to Improve Reproducibility and
Methodological Best Practices Adoption." _CRAN_.
<https://easystats.github.io/report/>.
  - Müller K, Wickham H (2023). _tibble: Simple Data Frames_. R package version
3.2.1, <https://CRAN.R-project.org/package=tibble>.
  - Patil I, Makowski D, Ben-Shachar M, Wiernik B, Bacher E, Lüdecke D (2022).
"datawizard: An R Package for Easy Data Preparation and Statistical
Transformations." _Journal of Open Source Software_, *7*(78), 4684.
doi:10.21105/joss.04684 <https://doi.org/10.21105/joss.04684>.
  - R Core Team (2023). _R: A Language and Environment for Statistical
Computing_. R Foundation for Statistical Computing, Vienna, Austria.
<https://www.R-project.org/>.
  - Robinson D, Hayes A, Couch S (2023). _broom: Convert Statistical Objects into
Tidy Tibbles_. R package version 1.0.5,
<https://CRAN.R-project.org/package=broom>.
  - Tang Y, Horikoshi M, Li W (2016). "ggfortify: Unified Interface to Visualize
Statistical Result of Popular R Packages." _The R Journal_, *8*(2), 474-485.
doi:10.32614/RJ-2016-060 <https://doi.org/10.32614/RJ-2016-060>,
<https://doi.org/10.32614/RJ-2016-060>. Horikoshi M, Tang Y (2018). _ggfortify:
Data Visualization Tools for Statistical Analysis Results_.
<https://CRAN.R-project.org/package=ggfortify>.
  - Wickham H (2016). _ggplot2: Elegant Graphics for Data Analysis_.
Springer-Verlag New York. ISBN 978-3-319-24277-4,
<https://ggplot2.tidyverse.org>.
  - Wickham H (2022). _stringr: Simple, Consistent Wrappers for Common String
Operations_. R package version 1.5.0,
<https://CRAN.R-project.org/package=stringr>.
  - Wickham H (2023). _forcats: Tools for Working with Categorical Variables
(Factors)_. R package version 1.0.0,
<https://CRAN.R-project.org/package=forcats>.
  - Wickham H (2023). _modelr: Modelling Functions that Work with the Pipe_. R
package version 0.1.11, <https://CRAN.R-project.org/package=modelr>.
  - Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G,
Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K,
Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K,
Yutani H (2019). "Welcome to the tidyverse." _Journal of Open Source Software_,
*4*(43), 1686. doi:10.21105/joss.01686 <https://doi.org/10.21105/joss.01686>.
  - Wickham H, François R, Henry L, Müller K, Vaughan D (2023). _dplyr: A Grammar
of Data Manipulation_. R package version 1.1.2,
<https://CRAN.R-project.org/package=dplyr>.
  - Wickham H, Henry L (2023). _purrr: Functional Programming Tools_. R package
version 1.0.1, <https://CRAN.R-project.org/package=purrr>.
  - Wickham H, Hester J, Bryan J (2023). _readr: Read Rectangular Text Data_. R
package version 2.1.4, <https://CRAN.R-project.org/package=readr>.
  - Wickham H, Vaughan D, Girlich M (2023). _tidyr: Tidy Messy Data_. R package
version 1.3.0, <https://CRAN.R-project.org/package=tidyr>.
  - Xie Y (2023). _knitr: A General-Purpose Package for Dynamic Report Generation
in R_. R package version 1.43, <https://yihui.org/knitr/>. Xie Y (2015).
_Dynamic Documents with R and knitr_, 2nd edition. Chapman and Hall/CRC, Boca
Raton, Florida. ISBN 978-1498716963, <https://yihui.org/knitr/>. Xie Y (2014).
"knitr: A Comprehensive Tool for Reproducible Research in R." In Stodden V,
Leisch F, Peng RD (eds.), _Implementing Reproducible Computational Research_.
Chapman and Hall/CRC. ISBN 978-1466561595.
report::report(fert_glmmTMB)
Random effect variances not available. Returned R2 does not account for random effects.
Warning in text == "" | text2 == "": longer object length is not a multiple of
shorter object length

Warning in text == "" | text2 == "": longer object length is not a multiple of
shorter object length
Random effect variances not available. Returned R2 does not account for random effects.
We fitted a linear model (estimated using ML and nlminb optimizer) to predict
YIELD with FERTILIZER (formula: YIELD ~ FERTILIZER). The model's explanatory
power related to the fixed effects alone (marginal R2) is 1.00. The model's
intercept, corresponding to FERTILIZER = 0, is at 51.93 (95% CI [29.18, 74.69],
p < .001). Within this model:

  - The intercept is statistically significant and positive (beta = 0.81, 95% CI
[0.66, 0.96], p < .001; Std. beta = 6.36e-10, 95% CI [-0.16, 0.16])

Standardized parameters were obtained by fitting the model on a standardized
version of the dataset. 95% Confidence Intervals (CIs) and p-values were
computed using a Wald z-distribution approximation. and We fitted a linear
model (estimated using ML and nlminb optimizer) to predict YIELD with
FERTILIZER (formula: YIELD ~ FERTILIZER). The model's explanatory power related
to the fixed effects alone (marginal R2) is 1.00. The model's intercept,
corresponding to FERTILIZER = 0, is at 16.99 (95% CI [10.96, 26.34]). Within
this model:

  - The intercept is statistically significant and positive (beta = 0.81, 95% CI
[0.66, 0.96], p < .001; Std. beta = 6.36e-10, 95% CI [-0.16, 0.16])

Standardized parameters were obtained by fitting the model on a standardized
version of the dataset. 95% Confidence Intervals (CIs) and p-values were
computed using a Wald z-distribution approximation.
report::report_text(fert_glmmTMB)
Random effect variances not available. Returned R2 does not account for random effects.
Warning in text == "" | text2 == "": longer object length is not a multiple of
shorter object length

Warning in text == "" | text2 == "": longer object length is not a multiple of
shorter object length
Random effect variances not available. Returned R2 does not account for random effects.
We fitted a linear model (estimated using ML and nlminb optimizer) to predict
YIELD with FERTILIZER (formula: YIELD ~ FERTILIZER). The model's explanatory
power related to the fixed effects alone (marginal R2) is 1.00. The model's
intercept, corresponding to FERTILIZER = 0, is at 51.93 (95% CI [29.18, 74.69],
p < .001). Within this model:

  - The intercept is statistically significant and positive (beta = 0.81, 95% CI
[0.66, 0.96], p < .001; Std. beta = 6.36e-10, 95% CI [-0.16, 0.16])

Standardized parameters were obtained by fitting the model on a standardized
version of the dataset. 95% Confidence Intervals (CIs) and p-values were
computed using a Wald z-distribution approximation. and We fitted a linear
model (estimated using ML and nlminb optimizer) to predict YIELD with
FERTILIZER (formula: YIELD ~ FERTILIZER). The model's explanatory power related
to the fixed effects alone (marginal R2) is 1.00. The model's intercept,
corresponding to FERTILIZER = 0, is at 16.99 (95% CI [10.96, 26.34]). Within
this model:

  - The intercept is statistically significant and positive (beta = 0.81, 95% CI
[0.66, 0.96], p < .001; Std. beta = 6.36e-10, 95% CI [-0.16, 0.16])

Standardized parameters were obtained by fitting the model on a standardized
version of the dataset. 95% Confidence Intervals (CIs) and p-values were
computed using a Wald z-distribution approximation.
report::report_model(fert_glmmTMB)
Random effect variances not available. Returned R2 does not account for random effects.
linear model (estimated using ML and nlminb optimizer) to predict YIELD with FERTILIZER (formula: YIELD ~ FERTILIZER)
report::report_parameters(fert_glmmTMB)
Random effect variances not available. Returned R2 does not account for random effects.
Warning in text == "" | text2 == "": longer object length is not a multiple of
shorter object length

Warning in text == "" | text2 == "": longer object length is not a multiple of
shorter object length
  - The intercept is statistically significant and positive (beta = 51.93, 95% CI [29.18, 74.69], p < .001; Std. beta = 6.36e-10, 95% CI [-0.16, 0.16])
  - The effect of FERTILIZER is statistically NA and positive (beta = 16.99, 95% CI [10.96, 26.34]; Std. beta = 0.96, 95% CI [0.79, 1.13])
  - The intercept is statistically significant and positive (beta = 0.81, 95% CI [0.66, 0.96], p < .001; Std. beta = 6.36e-10, 95% CI [-0.16, 0.16])
report::report_statistics(fert_glmmTMB)
Random effect variances not available. Returned R2 does not account for random effects.
Warning in text == "" | text2 == "": longer object length is not a multiple of
shorter object length

Warning in text == "" | text2 == "": longer object length is not a multiple of
shorter object length
beta = 51.93, 95% CI [29.18, 74.69], p < .001; Std. beta = 6.36e-10, 95% CI [-0.16, 0.16]
beta = 16.99, 95% CI [10.96, 26.34]; Std. beta = 0.96, 95% CI [0.79, 1.13]
beta = 0.81, 95% CI [0.66, 0.96], p < .001; Std. beta = 6.36e-10, 95% CI [-0.16, 0.16]
report::report_info(fert_glmmTMB)
Random effect variances not available. Returned R2 does not account for random effects.
Warning in text == "" | text2 == "": longer object length is not a multiple of
shorter object length

Warning in text == "" | text2 == "": longer object length is not a multiple of
shorter object length
Standardized parameters were obtained by fitting the model on a standardized version of the dataset. 95% Confidence Intervals (CIs) and p-values were computed using a Wald z-distribution approximation.
report::report_performance(fert_glmmTMB)
Random effect variances not available. Returned R2 does not account for random effects.
The model's explanatory power related to the fixed effects alone (marginal R2)
is 1.00
report::report_system()
Analyses were conducted using the R Statistical language (version 4.3.1; R Core
Team, 2023) on ArcoLinux
report::report(sessionInfo())
Analyses were conducted using the R Statistical language (version 4.3.1; R Core
Team, 2023) on ArcoLinux, using the packages effectsize (version 0.8.3;
Ben-Shachar MS et al., 2020), glmmTMB (version 1.1.7; Brooks ME et al., 2017),
effects (version 4.2.2; Fox J, Weisberg S, 2019), car (version 3.1.2; Fox J,
Weisberg S, 2019), carData (version 3.0.5; Fox J et al., 2022), lubridate
(version 1.9.2; Grolemund G, Wickham H, 2011), DHARMa (version 0.4.6; Hartig F,
2022), lindia (version 0.9; Lee Y, Ventura S, 2017), emmeans (version 1.8.7;
Lenth R, 2023), ggeffects (version 1.2.3; Lüdecke D, 2018), sjPlot (version
2.8.14; Lüdecke D, 2023), parameters (version 0.21.1; Lüdecke D et al., 2020),
performance (version 0.10.4; Lüdecke D et al., 2021), easystats (version 0.6.0;
Lüdecke D et al., 2022), see (version 0.7.4; Lüdecke D et al., 2021), insight
(version 0.19.3; Lüdecke D et al., 2019), bayestestR (version 0.13.0; Makowski
D et al., 2019), modelbased (version 0.8.6; Makowski D et al., 2020),
correlation (version 0.8.3; Makowski D et al., 2020), report (version 0.5.7;
Makowski D et al., 2023), tibble (version 3.2.1; Müller K, Wickham H, 2023),
datawizard (version 0.8.0; Patil I et al., 2022), broom (version 1.0.5;
Robinson D et al., 2023), ggfortify (version 0.4.16; Tang Y et al., 2016),
ggplot2 (version 3.4.2; Wickham H, 2016), stringr (version 1.5.0; Wickham H,
2022), forcats (version 1.0.0; Wickham H, 2023), modelr (version 0.1.11;
Wickham H, 2023), tidyverse (version 2.0.0; Wickham H et al., 2019), dplyr
(version 1.1.2; Wickham H et al., 2023), purrr (version 1.0.1; Wickham H, Henry
L, 2023), readr (version 2.1.4; Wickham H et al., 2023), tidyr (version 1.3.0;
Wickham H et al., 2023) and knitr (version 1.43; Xie Y, 2023).

References
----------
  - Ben-Shachar MS, Lüdecke D, Makowski D (2020). "effectsize: Estimation of
Effect Size Indices and Standardized Parameters." _Journal of Open Source
Software_, *5*(56), 2815. doi:10.21105/joss.02815
<https://doi.org/10.21105/joss.02815>, <https://doi.org/10.21105/joss.02815>.
  - Brooks ME, Kristensen K, van Benthem KJ, Magnusson A, Berg CW, Nielsen A,
Skaug HJ, Maechler M, Bolker BM (2017). "glmmTMB Balances Speed and Flexibility
Among Packages for Zero-inflated Generalized Linear Mixed Modeling." _The R
Journal_, *9*(2), 378-400. doi:10.32614/RJ-2017-066
<https://doi.org/10.32614/RJ-2017-066>.
  - Fox J, Weisberg S (2019). _An R Companion to Applied Regression_, 3rd
edition. Sage, Thousand Oaks CA.
<https://socialsciences.mcmaster.ca/jfox/Books/Companion/index.html>. Fox J,
Weisberg S (2018). "Visualizing Fit and Lack of Fit in Complex Regression
Models with Predictor Effect Plots and Partial Residuals." _Journal of
Statistical Software_, *87*(9), 1-27. doi:10.18637/jss.v087.i09
<https://doi.org/10.18637/jss.v087.i09>. Fox J (2003). "Effect Displays in R
for Generalised Linear Models." _Journal of Statistical Software_, *8*(15),
1-27. doi:10.18637/jss.v008.i15 <https://doi.org/10.18637/jss.v008.i15>. Fox J,
Hong J (2009). "Effect Displays in R for Multinomial and Proportional-Odds
Logit Models: Extensions to the effects Package." _Journal of Statistical
Software_, *32*(1), 1-24. doi:10.18637/jss.v032.i01
<https://doi.org/10.18637/jss.v032.i01>.
  - Fox J, Weisberg S (2019). _An R Companion to Applied Regression_, Third
edition. Sage, Thousand Oaks CA.
<https://socialsciences.mcmaster.ca/jfox/Books/Companion/>.
  - Fox J, Weisberg S, Price B (2022). _carData: Companion to Applied Regression
Data Sets_. R package version 3.0-5,
<https://CRAN.R-project.org/package=carData>.
  - Grolemund G, Wickham H (2011). "Dates and Times Made Easy with lubridate."
_Journal of Statistical Software_, *40*(3), 1-25.
<https://www.jstatsoft.org/v40/i03/>.
  - Hartig F (2022). _DHARMa: Residual Diagnostics for Hierarchical (Multi-Level
/ Mixed) Regression Models_. R package version 0.4.6,
<https://CRAN.R-project.org/package=DHARMa>.
  - Lee Y, Ventura S (2017). _lindia: Automated Linear Regression Diagnostic_. R
package version 0.9, <https://CRAN.R-project.org/package=lindia>.
  - Lenth R (2023). _emmeans: Estimated Marginal Means, aka Least-Squares Means_.
R package version 1.8.7, <https://CRAN.R-project.org/package=emmeans>.
  - Lüdecke D (2018). "ggeffects: Tidy Data Frames of Marginal Effects from
Regression Models." _Journal of Open Source Software_, *3*(26), 772.
doi:10.21105/joss.00772 <https://doi.org/10.21105/joss.00772>.
  - Lüdecke D (2023). _sjPlot: Data Visualization for Statistics in Social
Science_. R package version 2.8.14,
<https://CRAN.R-project.org/package=sjPlot>.
  - Lüdecke D, Ben-Shachar M, Patil I, Makowski D (2020). "Extracting, Computing
and Exploring the Parameters of Statistical Models using R." _Journal of Open
Source Software_, *5*(53), 2445. doi:10.21105/joss.02445
<https://doi.org/10.21105/joss.02445>.
  - Lüdecke D, Ben-Shachar M, Patil I, Waggoner P, Makowski D (2021).
"performance: An R Package for Assessment, Comparison and Testing of
Statistical Models." _Journal of Open Source Software_, *6*(60), 3139.
doi:10.21105/joss.03139 <https://doi.org/10.21105/joss.03139>.
  - Lüdecke D, Ben-Shachar M, Patil I, Wiernik B, Makowski D (2022). "easystats:
Framework for Easy Statistical Modeling, Visualization, and Reporting." _CRAN_.
R package, <https://easystats.github.io/easystats/>.
  - Lüdecke D, Patil I, Ben-Shachar M, Wiernik B, Waggoner P, Makowski D (2021).
"see: An R Package for Visualizing Statistical Models." _Journal of Open Source
Software_, *6*(64), 3393. doi:10.21105/joss.03393
<https://doi.org/10.21105/joss.03393>.
  - Lüdecke D, Waggoner P, Makowski D (2019). "insight: A Unified Interface to
Access Information from Model Objects in R." _Journal of Open Source Software_,
*4*(38), 1412. doi:10.21105/joss.01412 <https://doi.org/10.21105/joss.01412>.
  - Makowski D, Ben-Shachar M, Lüdecke D (2019). "bayestestR: Describing Effects
and their Uncertainty, Existence and Significance within the Bayesian
Framework." _Journal of Open Source Software_, *4*(40), 1541.
doi:10.21105/joss.01541 <https://doi.org/10.21105/joss.01541>,
<https://joss.theoj.org/papers/10.21105/joss.01541>.
  - Makowski D, Ben-Shachar M, Patil I, Lüdecke D (2020). "Estimation of
Model-Based Predictions, Contrasts and Means." _CRAN_.
<https://github.com/easystats/modelbased>.
  - Makowski D, Ben-Shachar M, Patil I, Lüdecke D (2020). "Methods and Algorithms
for Correlation Analysis in R." _Journal of Open Source Software_, *5*(51),
2306. doi:10.21105/joss.02306 <https://doi.org/10.21105/joss.02306>,
<https://joss.theoj.org/papers/10.21105/joss.02306>.
  - Makowski D, Lüdecke D, Patil I, Thériault R, Ben-Shachar M, Wiernik B (2023).
"Automated Results Reporting as a Practical Tool to Improve Reproducibility and
Methodological Best Practices Adoption." _CRAN_.
<https://easystats.github.io/report/>.
  - Müller K, Wickham H (2023). _tibble: Simple Data Frames_. R package version
3.2.1, <https://CRAN.R-project.org/package=tibble>.
  - Patil I, Makowski D, Ben-Shachar M, Wiernik B, Bacher E, Lüdecke D (2022).
"datawizard: An R Package for Easy Data Preparation and Statistical
Transformations." _Journal of Open Source Software_, *7*(78), 4684.
doi:10.21105/joss.04684 <https://doi.org/10.21105/joss.04684>.
  - R Core Team (2023). _R: A Language and Environment for Statistical
Computing_. R Foundation for Statistical Computing, Vienna, Austria.
<https://www.R-project.org/>.
  - Robinson D, Hayes A, Couch S (2023). _broom: Convert Statistical Objects into
Tidy Tibbles_. R package version 1.0.5,
<https://CRAN.R-project.org/package=broom>.
  - Tang Y, Horikoshi M, Li W (2016). "ggfortify: Unified Interface to Visualize
Statistical Result of Popular R Packages." _The R Journal_, *8*(2), 474-485.
doi:10.32614/RJ-2016-060 <https://doi.org/10.32614/RJ-2016-060>,
<https://doi.org/10.32614/RJ-2016-060>. Horikoshi M, Tang Y (2018). _ggfortify:
Data Visualization Tools for Statistical Analysis Results_.
<https://CRAN.R-project.org/package=ggfortify>.
  - Wickham H (2016). _ggplot2: Elegant Graphics for Data Analysis_.
Springer-Verlag New York. ISBN 978-3-319-24277-4,
<https://ggplot2.tidyverse.org>.
  - Wickham H (2022). _stringr: Simple, Consistent Wrappers for Common String
Operations_. R package version 1.5.0,
<https://CRAN.R-project.org/package=stringr>.
  - Wickham H (2023). _forcats: Tools for Working with Categorical Variables
(Factors)_. R package version 1.0.0,
<https://CRAN.R-project.org/package=forcats>.
  - Wickham H (2023). _modelr: Modelling Functions that Work with the Pipe_. R
package version 0.1.11, <https://CRAN.R-project.org/package=modelr>.
  - Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R, Grolemund G,
Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller E, Bache SM, Müller K,
Ooms J, Robinson D, Seidel DP, Spinu V, Takahashi K, Vaughan D, Wilke C, Woo K,
Yutani H (2019). "Welcome to the tidyverse." _Journal of Open Source Software_,
*4*(43), 1686. doi:10.21105/joss.01686 <https://doi.org/10.21105/joss.01686>.
  - Wickham H, François R, Henry L, Müller K, Vaughan D (2023). _dplyr: A Grammar
of Data Manipulation_. R package version 1.1.2,
<https://CRAN.R-project.org/package=dplyr>.
  - Wickham H, Henry L (2023). _purrr: Functional Programming Tools_. R package
version 1.0.1, <https://CRAN.R-project.org/package=purrr>.
  - Wickham H, Hester J, Bryan J (2023). _readr: Read Rectangular Text Data_. R
package version 2.1.4, <https://CRAN.R-project.org/package=readr>.
  - Wickham H, Vaughan D, Girlich M (2023). _tidyr: Tidy Messy Data_. R package
version 1.3.0, <https://CRAN.R-project.org/package=tidyr>.
  - Xie Y (2023). _knitr: A General-Purpose Package for Dynamic Report Generation
in R_. R package version 1.43, <https://yihui.org/knitr/>. Xie Y (2015).
_Dynamic Documents with R and knitr_, 2nd edition. Chapman and Hall/CRC, Boca
Raton, Florida. ISBN 978-1498716963, <https://yihui.org/knitr/>. Xie Y (2014).
"knitr: A Comprehensive Tool for Reproducible Research in R." In Stodden V,
Leisch F, Peng RD (eds.), _Implementing Reproducible Computational Research_.
Chapman and Hall/CRC. ISBN 978-1466561595.

References

Fowler, J., L. Cohen, and P. Jarvis. 1998. Practical Statistics for Field Biology. England: John Wiley & Sons.